% Define file paths
filenames = {
    "C:\Users\86180\Desktop\beamA\data\filter\483 filter spectrum.txt",
    "C:\Users\86180\Desktop\beamA\data\filter\704 filter spectrum.txt",
    "C:\Users\86180\Desktop\beamA\data\filter\light source spectrum.txt",
    "C:\Users\86180\Desktop\beamA\data\filter\beam spliter.txt",
    "C:\Users\86180\Desktop\beamA\data\filter\light source spectrum for beam spliter.txt"
};

% Initialize structure to store spectral data
spectraData = struct();

% Load data from files
for i = 1:length(filenames)
    fileID = fopen(filenames{i}, 'r');
    if fileID == -1
        fprintf('Cannot open file: %s\n', filenames{i});  
        continue;
    end

    while ~feof(fileID)
        tline = fgetl(fileID);
        if contains(tline, '>>>>>Begin Processed Spectral Data<<<<<')
            break;
        end
    end

    data = textscan(fileID, '%f %f', 'Delimiter', '\t');
    spectraData(i).wavelength = data{1};
    spectraData(i).intensity = data{2};
    
    fclose(fileID);

    if isempty(spectraData(i).wavelength)
        fprintf('No data: %s\n', filenames{i}); 
    else
        fprintf('Data read successfully: %s with %d data points.\n', filenames{i}, length(spectraData(i).wavelength)); 
    end
end

% Analyze and plot spectra data for filter 483 and 704 with light source spectrum
% Analyze and plot spectra data for filter 483 and 704 with light source spectrum
for i = 1:2  % Only process the first two spectra for 483 and 704 filters
    figure(i);
    hold on;
    
    % Correct way to set filterLabel based on index
    if i == 1
        filterLabel = '483 filter';
    else
        filterLabel = '704 filter';
    end
    
    plot(spectraData(i).wavelength, spectraData(i).intensity, 'b-', 'DisplayName', filterLabel);
    plot(spectraData(3).wavelength, spectraData(3).intensity, 'r-', 'DisplayName', 'Light source');
    xlabel('Wavelength (nm)');
    ylabel('Intensity');
    title(sprintf('%s Spectrum with Light Source', filterLabel));
    legend show;
    
    % Find significant peak and FWHM
    [peakWavelength, fwhmWidth] = robustPeakFindAndFWHM(spectraData(i).wavelength, spectraData(i).intensity);
    % Calculate light intensity ratio at the peak
    if ~isnan(peakWavelength)
        lightSourceIntensity = interp1(spectraData(3).wavelength, spectraData(3).intensity, peakWavelength, 'linear', 'extrap');
        filterIntensity = interp1(spectraData(i).wavelength, spectraData(i).intensity, peakWavelength, 'linear', 'extrap');
        intensityRatio = lightSourceIntensity / filterIntensity;
        fprintf('Figure %d - %s: Peak at %.2f nm with FWHM %.2f nm. Intensity Ratio (Light Source/Filter): %.2f\n', i, filterLabel, peakWavelength, fwhmWidth, intensityRatio);
    end

    hold off;
end


% Function to find significant peaks and calculate FWHM
function [peakWavelength, fwhmWidth] = robustPeakFindAndFWHM(x, y)
    [pks, locs, widths, proms] = findpeaks(y, x, 'MinPeakProminence', max(y) * 0.1);
    if isempty(pks)
        peakWavelength = NaN;
        fwhmWidth = NaN;
        return;
    end
    [~, idxMaxProm] = max(proms);
    peakWavelength = locs(idxMaxProm);
    peakIntensity = pks(idxMaxProm);
    halfMax = peakIntensity / 2;
    leftIdx = find(y(1:idxMaxProm) <= halfMax, 1, 'last');
    rightIdx = find(y(idxMaxProm:end) <= halfMax, 1, 'first') + idxMaxProm - 1;
    if isempty(leftIdx) || isempty(rightIdx)
        fwhmWidth = NaN;
    else
        x1 = interp1(y([leftIdx leftIdx+1]), x([leftIdx leftIdx+1]), halfMax, 'linear', 'extrap');
        x2 = interp1(y([rightIdx rightIdx+1]), x([rightIdx rightIdx+1]), halfMax, 'linear', 'extrap');
        fwhmWidth = x2 - x1;
    end
end








